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A mathematical model describing a steady pH-gradient in the solution of ampholytes in water has 
been studied with the use of analytical, asymptotic, and numerical methods. We show that at the 
large values of an electric current a concentration distribution takes the form of a piecewise constant 
function that is drastically different from a classical Gaussian form. The correspondent pH-gradient 
takes a stepwise form, instead of being a linear function. A discovered anomalous pH-gradient can 
crucially affect the understanding of an isoelectric focusing process. 
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I. INTRODUCTION 



The isoelectric focusing (IEF) is a classical method of 
a mixture separation that is widely used in medicine and 
biology. There is an enormous number of publications 
devoted to its theory and applications. We mention here 
only fundamental papers [3,11, 1 0, Hi, the monographs 
considering various theoretical and practical aspects of 
the IEF_H0,J1 , and papers [1 El, El El El, El, El, 
E1E3 1 18L 1 19L [20( that contain mathematical models and 
computer simulations that are close to the topic of our 
paper. There are also a number of papers that study the 
relatedproblems in space biotechnology [U, [H, [H, l24l . 

EEUE 

An important part of the IEF is the creation of stable 
pH-gradients that are used for the setting up of a focusing 
process. Carrier ampholytes have been in a wide prac- 
tical use for the creation of pH-gradients starting from 
pioneering papers 

Q, lH lH (see also @,E3)- 

The use 

of immobilized pH-gradients for the same purpose has 
been proposed and introduced by 0, EH , while the use 
of borate-polyol systems for the creation of pH-gradients 
is considered in (28[. The mathematical models describ- 
ing the formation of stable pH-gradients with the use 
of carrier ampholytes are given in [H, El, E3, El, HU- 
In particular, [ll, [2I] introduce the notion of mixtures 
with the infinite number of components (the infinite- 
component mixtures) where the authors formally replace 
the discrete number k in the concentration Ck onto the 
continuous parameter s; in other words the discrete set 
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Cfc, k = 1, ...,n is replaced by the distribution func- 
tion c(s). The first and simplest mathematical models 
describing the creation of linear (or almost linear) pH- 
gradients with the use of carrier ampholytes were given 
in 0,0] • More advanced approaches that employ the gov- 
erning differential equations of continuous media, mathe- 
matical modeling, and computations are given in [l5l.[la |. 
These papers show that an ampholyte concentration dis- 
tribution should be close to Gaussian's one and a pH- 
gradicnt should be linear. Similar results with the use of 
models for infinite-component mixtures were obtained in 

The statement that a concentration distribution pos- 
sesses a Gaussian shape is currently accepted as a key 
classical result. Moreover, it is widely assumed that the 
increase of a potential difference (or an electric current) 
stabilizes a linear pH-gradient and intensifies a focusing 
process. On the contrary, the main results of our paper 
show that the increase of an electric current (starting 
from some threshold value) causes an anomalous steady 
stepwise pH-gradient and a piecewise-constant concen- 
tration distribution. In addition we present explicit ana- 
lytical expressions for a pH distribution and for an elec- 
tric conductivity as the functions of ampholyte concen- 
trations. The obtained asymptotic expressions for the 
concentrations of ampholytes allow us to carry out the 
calculations of stepwise pH-gradients. 



II. MATHEMATICAL MODEL 

Well-known governing equations describing IEF in the 
aqueous solution of ampholyts in dimensionless variables 
are El El, Hi: 



dc k 
dt 



dik 
dx 



0, k = l, 



(1) 



2 



ik = ~ £ ^~q^ +lkek(ip)ckE, (2) 



In this case the concentrations of ions H + and OH~ 
reflated as 



are 



\- / d{e k {ip)c k/ 
J = l^lkAs + a k (ip)c k E 

k=l 



e k (tp)ck + 2K W sinh ip = 0, 



fc=i 



where 



1 -f 
* 7h 



(4) 



(5) 



lgh = -lgK w -ipige. (6) 



Here 0^(2;, t) are molar ampholyte concentrations, ip(x, t) 
is a function, linearly connected to pH of a solution, 
E(x, t) is the strength of an electric field, i k are the densi- 
ties of molar ampholyte fluxes, j = const is the density of 
an electric current (it can also depend on time j = j(t)), 
his a hydrogen ion concentration, e is a diffusion parame- 
ter, 7^ > are characteristic mobilities of components in 
an electric field, ej k > are diffusion coefficients, e k {ip), 
o~ k (ip) are the molar charges and the molar conductivities 
of ampholytes, "f k e k (ip) are the electrophoretic mobilities 
of ampholytes, "/h, Ioh are the mobilities of hydrogen 
H + and hydroxyl OEP ions, 70 is the effective mobility 
of water ions, K w is the autodissociation constant for 
water, ipo is the value of ip when water conductivity is 
minimal. 

The system of equations ©-© describes the mass 
transport under the action of an electric field in a chem- 
ically active medium. Chemical reactions in the medium 
are assumed to be very fast (being completed almost in- 
stantly). The functions e k (ip), cr k {ip) are given; in each 
particular case they can be defined from the chemical ki- 
netic equations (see below). The algebraic equation © 
represents the condition of mixture electroneutrality; it 
allows us to find ip. In fact, this equation describes the 
instant control of medium properties (electrophoretic mo- 
bilities and molar conductivities) by the function ip (that 
is linked to the concentration of hydrogen ions). Such 
situation is typical only for a weak aqueous solution; in 
more general cases the values e k , a k can depend on c k . 
The term 2i^„,7 (E — edip/dx) cosh(-0 — ip ) m © de- 
scribes the contribution of water ions into mixture con- 
ductivity, while the term 2K W sinh ip in © corresponds 
to the contribution of these ions into a mixture molar 
charge. These terms should be included into the gov- 
erning equations if the autodissociation reaction is taken 
into account: 



H 2 O^H + + OIT 

k' 



h=[H+], [OH-] = ^ = 



K 2 
h 



where is an autodissociation constant for water, 
square brackets [ • ] denote the values of concentrations, 
as it is accepted in chemistry. 

Let us note that the maximal values of concentrations 
Cfe and the values 70 , ipo , j k , ip k , k = I , . . . , n are all of 
the order O(l), while the parameters e and K w are small. 
Notice, that we denote an autodissociation constant for 
water as , not K w ; hence K w represents a dissociation 
constant for water, its value in dimensional variables is 
K w = 10 _7 mol/l. In addition, instead of a convention- 
ally used function pH (measure of the acidity or alkalinity 
of a solution) we use function ip (linearly connected with 
pH) which is better adapted to our mathematical model 
(see ©) H. 

In order to make the system (U}-© complete, we have 
to prescribe the functions e k (ip), cr k (ip) for ampholytes. 
Let us consider the solution of n amphoteric substances 
(ampholytes) that exhibit both acid and based properties 
(a typical representative is an amino acid). The dissoci- 
ation reactions are (i = 1, . . . , n) [H, I20I [29j: 



H+R ^ R? 



IT 



A i 

R ,- ^ R , 



FT 



(7) 



Here R° is a zwitterion (for instance an amino acid 
residue), Ai and Bi are dissociation constants for acid 
(R^) an d based (H + Ri) groups correspondingly. Notice, 
that even for polyampholytes with supplementary side 
groups a reaction scheme can be described by ©, as- 
suming the pH of a solution being close to the isoelectric 
point of a substance Q. The analytical concentrations 
q, molar charges ei{ip), and molar conductivities <7i(ip) 
for © are: 



Ci 



[H+R,] + [R°] + [Rr], 



1 dtpijiP) I d 2 tpi(i>) 

e AW) = 7TT , <W) " 



y>i(ip) dip 



ipi(ip) dip 2 



tfiiiip) = Si + cosh(ip - ipi), 



(8) 



1, MBi 
— In r— 

2 K 2 



1 (Ek 

2 I Ai 



Here ip k is an isoelectric point (the value of ip when the 
molar charge e k is zero, i.e. e k (ip k ) = 0), tfi and Si 
are auxiliary notations. Notice, that it is convenient not 
to write explicit expressions for e k (ip) and o~ k (ip) but to 
present them as the derivatives of functions ifi(ip). 
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III. FORMULATION OF STEADY 
PH-GRADIENT PROBLEM 

Let us take dc k / dt = and introduce non-leak bound- 
ary conditions for the ampholytes at the walls of an elec- 
trophorctic chamber (0 ^ x ^ 1): 

i fc (0) = 0, * fc (l) = 0, k = l,...,n. (9) 

The total molar concentration of each ampholyte is: 



c k (x)dx = M fe , k = 1, 



(10) 



The values 70, V'o, Ik, ipk, M k , k = 1, . . . , n, K w , e are 
given. For the sake of simplicity we consider only the 
IEF case with a given electric current density j = const. 
The problem with a given constant voltage can be solved 
similarly however it is more tedious. The integration of 
equations ([I]) with boundary conditions (|9|) yields i k = 
that can be written as: 



dck 



dx 



e k {i>)c k E = 0, k = 1 
The remaining equations keep their form: 

_d(e k (ijj)c k ) 



(11) 



3 = J2^ k 

fc=i 



— s- 



dx 



a k (ip)c k E 



dtp 



2K w j Q [ E - e— ) cosh('0 - ^0), 



e k {ip)c k + 2K W smhip = 0. 



(12) 



(13) 



k=l 



The system of equations (fTT j) - (fl~3"| with additional con- 
ditions (fTU|> allow to find the concentrations c k (x), k = 
1, . . . , n, the strengths of the electric field E(x), and the 
function tp(x) (i.e. pH(x)). Let us emphasize here that 
despite of the presence of derivative dtp/dx we do not 
need to introduce any supplementary condition for "0, 
since the function ij) has been determined by the alge- 
braic equation (fl"3"l) . 



IV. 



THE TRANSFORMATION OF THE 
PROBLEM 



For the problem (fTTTj) — (jT5|) we have found the change 
of dependent variables that allows us to obtain explicit 
analytical expressions for functions ip{ x ) an( i E(x). We 
introduce a new 'concentration' a k (x) and the 'effective' 
strength of the electric field Eq(x) as: 



c k (x) = a k (x)ip k (i/>(x)), 



E (x) = E(x) -e 



dip(x) 
dx 



(14) 



In this case one can easily show that il>( x ) an d E[x) are 
defined by the explicit expressions, provided that func- 
tions a k (x) are known: 



2K W + £ a k (x)e* 



ip(x) = - In 



k=l 



2K W + a k (x)e-4* 

k=l 



E(x)=E (x)+e^±, E (x) = ^-, (15) 



dx 



a{x) 



where 



- / d 2 \ 

o-(^) = 22lkak{x) yPkii')-^ ^ a - ( Ph{i>)j 

+ 2K w j cosh(^(a;) - Vo)- 



(16) 



In order to find a k (x) (taking into account the conditions 
(flC)]) and the changes of variables (fl4|) ') we have the fol- 
lowing system of 2n differential equations with ij) and a 
defined by {HJ and (16]) 



dx k ^ ^ cr(x) dip 



— In f k (ip) 



(17) 



— m k (x) = a k (p k (ip(x)), k = l,...,n. (18) 
dx 



£0 



J 



m k (x) = \ Cfe(s) ds 



with the boundary conditions: 

m fe (0) = 0, m k (l) = M k , k = l, 



,n. (19) 



The additional variables m k have the meaning of the total 
number of moles of an ampholyte in the interval [0, x]. 



V. THE MAIN ASYMPTOTIC TERM AT K w -> 

The nonlinear boundary problem (fT7|) — (|T5]) can be 
solved numerically with the use of a shooting method, 
by the reduction of (fT7|) - (fT9j) to Cauchy's problem with 
unknown initial data a k (0). To find a k (0) one should 
solve n nonlinear transcendent equations. It can be per- 
formed for example by the Newton-Raphson method that 
requires good initial approximations. An initial approxi- 
mation for small £0 and K w can be obtained with the use 
of asymptotic formulas. The analysis of equations (JTTJ) 
show that their right side parts at Eq — > are propor- 
tional to dip k (ip)/dip ~ sinh(?/' — ip k ) (see (0). It means 
that ip(x) — ► ip k for k = 1, . . . , n. It is possible only if 
functions a k (x), k = l,...,n, ip(x), and E(x) are close 
to stepwise ones. 
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Let us introduce the set of intervals (xk-i, x k ) G [0, 1] 
that are not overlapping. We assume that 

Ofc(x) = a°, xe(x k -i,x k ) k = l,...,n, 



a k (x)=0, x $ (x k -i,x k ), k = l,...,n, (20) 

where are some constants. Then a particular conclu- 
sion from (|15[) at K w — > is: 

ip(x) = tp k , x€(x k -i,x k ), k = l, ...,n. (21) 

Next we introduce the functions: 

1 1 d^fcMaO) 



$fc(z) 



cr(x) ip k (ifj(x)) dijj 



X 

F k (x) = — <P k (s)ds, k = l,...,n. 
£q J 



For an interval (xi-i,Xi) we get: 



(22) 



7^ 



ft ' 



(23) 



eoF fc (a;) = ^ $^(x s - z s _i) + ^\x - xt-x). 



The integrating of (JT7J) and taking into account notations 
([22]) yield: 



a fc (x) = a fe (0)e F *H 



(24) 



Using (El), and ([21]), the equations a k (x) = a° k for 

x G (x fe _i,a;fc) give: 



(0) = a°exp (~ y E*t\x.-x a - 1 )) . (25) 



s=l 



The values a® can be easily found from the boundary 
condition m k (l) = M k with the help of ([20]). (|2Tj) : 



«°ft 



(4 + l)(x fe - x k -x) ' 



k = l,...,n. (26) 



Expressions ([26]) . (f25j) . ([23]) represent the required ini- 
tial approximation for a k (0), provided that the lengths 
of intervals (x k , x k _\), k = 1, . . . ,n are given. The main 
asymptotic term at K w — > is given by 



VI. NUMERICAL RESULTS 

The main difficulty in the use of formulas (f2"5|) 
is the fact that the lengths of the intervals {x k -i,x k ), 



k = 1 , . . . , n are unknown, they can be determined by the 
next asymptotic terms. However our numerical experi- 
ments have shown that (at least for n > 5) it is sufficient 
to take their lengths equal: (x k — x k ^\) = 1/n. We have 
shown that asymptotic formulas and numerical results 
are in good agreement beginning from £o ~ 2 • 1CP 4 and 
K w sa 10~ 4 . In particular, the numerical results obtained 
by a shooting method for the values of parameters n = 8, 
4 = 99, Vfc = 3.85 - 0.6fc, M k = 160, e = 5.25 • 10~ 6 , 
7 fc = 0.15, = -0.274, 70 = 7.077, K w = 10" 5 are 
given in Table U where for the sake of convenience we 
present the values z k = — lOeohiafc(O). The values z£ 



k 


2 


3 


4 


5 


6 


7 


8 


Zk 


0.0042 


0.0138 


0.0325 


0.0662 


0.1256 


0.2264 


0.3894 


Zfc 


0.0039 


0.0137 


0.0324 


0.0664 


0.1262 


0.2276 


0.3914 



TABLE I: The comparison between the numerical and asymp- 
totic values 



in the second raw correspond to their asymptotic values. 
One should notice that the formulas (f2"S")) . (|2"6")) do not 
work well for ai(0), so for the obtaining of its asymptotic 
value one has to use an additional procedure that are 
similar to Sect. II VI for a k (l). Fig. [T]shows the concentra- 
tion distributions a k (x) and function ip( x ) for t w ° values 
e = 5.25 • 10~ 5 (left) and e = 5.25 • 10~ 6 (right). It 
is clearly seen that a k (x) are located successively on the 
x-axis in the interval (x k _i, x k ). The concentration max- 
ima correspond to the isoelectric points i/j k . In a classical 
case when eg is not too small, function ip(x) is monotoni- 
cally decreasing and concentration distributions are close 
to Gaussian ones (the left frame). The decreasing of eq 
(or the increasing of current density) causes the 'almost' 
stepwise distribution of ip(x) i.e. ip{x) ~ i\) k in the in- 
tervals (x k -i,x k ) and 'almost' constant concentrations 



a k (x) 



(the right frame). 
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